Introduction

General remarks

The R pacakge enerscape computes the energy landscape based on models predicting the energy costs of transport. Two models are currently implemented, the ARC model for terrestrial, legged animals from (Pontzer 2016) and the model for cycling expenditures from (Prampero et al. 1979). Adding other models is relatively easy and can be done by writing a custom function in the enerscape_internals.R script. Currently, there is no implementation for user-specified models, but I plan to add this in the future.

An extensive derivation of the ARC model can be found in Pontzer (2016). Here, it is sufficient to say that the ARC model accounts for the energy spent for cross-bridge cycling, generating tension, as well as associated activation-relaxation processes that restore cellular ion gradients - mostly Ca\(^{2+}\) necessary for cross-bridge cycling. Notably, predictions of the ARC model fit very well observed data. In enerscape, the dataset from Pontzer (2016) can be access as the variable pontzer:

library(tidyverse)
library(raster)
library(sf)
library(enerscape)

#' Function to plot the areas
#' 
#' @param x is a raster
#' @param poly is a shapefile polygon
plot_area <- function(x, poly, mask = NULL, col = NULL, void = FALSE) {
  if (is.null(col)) {
    plot(x, axes = !void, box = !void)
  } else {
    plot(x, col = col, axes = !void, box = !void)
  }
  if (!is.null(mask)) {
    plot(mask, add = T, legend = FALSE,
         col = adjustcolor("white", alpha.f = 0.5), 
         axes = !void, 
         box = !void)
  }
  plot(poly$geometry, add = TRUE)
}

#' ARC model predictions
#'
#' @param mass body mass of the animal (kg)
#' @param slope incline of the terrain (degree)
ARC <- function(mass, slope) {
  E_ar <- 8 * mass ^ (-0.34)
  E_mec <- 100 * (1 + sin((2 * slope - 74) / 180 * pi)) * mass ^ (-0.12)
  E <- (E_ar + E_mec)
  return(E)
}

pontzer %>% 
  as_tibble %>% 
  mutate(ARC = ARC(Mass, Incline)) %>% 
  ggplot() +
  geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
  geom_point(aes(Cost.of.Transport, ARC, col = Incline)) +
  scale_color_gradient(low = "steelblue", high = "tomato") +
  scale_x_log10() +
  scale_y_log10() +
  xlab(expression("E"["COT"]*" (J kg"^"-1"*"m"^"-1"*")")) +
  ylab("ARC predicitons") +
  theme_bw()

Computing energy landscapes

To compute energy landscapes using enerscape, we need an elevation rasterLayer, normally called a digital elevation model (DEM). The raster should have equal resolution on both horizontal and vertical axes and a planar coordinate reference systems in meters. The following is a raster that have these characteristics:

dem <- raster("../data/dem-vignette.tif")
res(dem)
## [1] 100 100
crs(dem)
## CRS arguments:
##  +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs

We can then compute the energy landscape of the area for a general animal of 50 kg using the enerscape() function:

en <- enerscape(dem, 50)
##  - Raster cells are assumed to have same horizontal and vertical resolution and with planar coordinate reference system (e.g. UTM)
##   | Calculating slope
##   | Calculating work
##   | Calculating conductance (1 / work)
##  - Do not use slope with negative values for distance calculations

The output is an enerscape object, which is a list with elements:

  • The neighbors used to compute the transition matrix among cells, either 4 (in chess, rook’s move of 1), 8 (king’s move), or 16 (king’s + knigth’s move).
  • The body mass of the animal (kg).
  • A rasterStack with DEM, Slope, Work (cost of travel), and Conductance (1 / Work).
  • A conductance transition matrix that summarize the conductance of moving between neighboring cells.

The function en_extrapolation() can be used to highlight where ARC predictions extrapolate from the testing data. A rasterLayer masking extrapolations for the incline is returned if available:

extr <- en_extrapolation(en, plot = FALSE)
extr
## $`Mass extrapolated`
## [1] FALSE
## 
## $`Slope extrapolated`
## [1] TRUE
## 
## $`Slope extrapolation`
## class      : RasterLayer 
## dimensions : 367, 509, 186803  (nrow, ncol, ncell)
## resolution : 100, 100  (x, y)
## extent     : 852140, 903040, 4659980, 4696680  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : Slope 
## values     : 1, 1  (min, max)
plot(dem)
plot(extr$`Slope extrapolation`, 
     add = TRUE, 
     col = adjustcolor("blue", alpha.f = 0.3),
     legend = FALSE)
arrows(888523, 4674275, 878791, 4671549)

The algorithm to compute least-cost paths from the R package gdistance is implemented in the en_lcp() functio0n (Etten 2017). en_lcp() takes as input two points for which the least-cost path need to be evaluated. Alternatively, random points are generated if simulate_random_points = TRUE:

lcp <- enerscape::en_lcp(en, simulate_random_points = TRUE)
## 1 of 10 random points ...
## 2 of 10 random points ...
## 3 of 10 random points ...
## 4 of 10 random points ...
## 5 of 10 random points ...
## 6 of 10 random points ...
## 7 of 10 random points ...
## 8 of 10 random points ...
## 9 of 10 random points ...
## 10 of 10 random points ...

The random walk algorithm from gdistance is also implemented as en_passage():

walk <- enerscape::en_passage(en,
                              simulate_random_points = TRUE,
                              theta = 4)
## 1 of 10random points ...
## 2 of 10random points ...
## 3 of 10random points ...
## 4 of 10random points ...
## 5 of 10random points ...
## 6 of 10random points ...
## 7 of 10random points ...
## 8 of 10random points ...
## 9 of 10random points ...
## 10 of 10random points ...

walk
## $Origins
##            x       y
##  [1,] 861890 4687030
##  [2,] 868490 4690430
##  [3,] 852890 4682130
##  [4,] 894290 4665930
##  [5,] 855990 4678530
##  [6,] 866890 4686230
##  [7,] 900490 4692930
##  [8,] 893790 4681730
##  [9,] 874190 4666630
## [10,] 857690 4681530
## 
## $Destinations
##            x       y
##  [1,] 902690 4690630
##  [2,] 852490 4695830
##  [3,] 863090 4677930
##  [4,] 864190 4674430
##  [5,] 895690 4679530
##  [6,] 869790 4667330
##  [7,] 889990 4686530
##  [8,] 902590 4666030
##  [9,] 855890 4679330
## [10,] 863790 4662530
## 
## $Passages
## class      : RasterStack 
## dimensions : 367, 509, 186803, 10  (nrow, ncol, ncell, nlayers)
## resolution : 100, 100  (x, y)
## extent     : 852140, 903040, 4659980, 4696680  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs 
## names      : layer.1, layer.2, layer.3, layer.4, layer.5, layer.6, layer.7, layer.8, layer.9, layer.10 
## min values :       0,       0,       0,       0,       0,       0,       0,       0,       0,        0 
## max values :       0,       0,       0,       0,       0,       0,       0,       0,       0,        0 
## 
## 
## $`Cumulative net passage`
## class      : RasterLayer 
## dimensions : 367, 509, 186803  (nrow, ncol, ncell)
## resolution : 100, 100  (x, y)
## extent     : 852140, 903040, 4659980, 4696680  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : NA, NA  (min, max)

As this may be quite slow, I implemented another way of computing the overall connectivity between two points - or for the whole landscape - using Circuitscape and Omniscape. The initialization files for these to algorithms can be written using circuitscape_skeleton() and omniscape_skeleton(), which create circuitscape.ini and omniscape.ini to be launched from within a Julia terminal. See the example [Landscape connectivity for the Marsican bear (Ursus arctos marsicanus) in the Sirente-Velino Regional Park: Circuitscape and Omniscape implementations.] for details.

Landscape connectivity for the Marsican bear (Ursus arctos marsicanus) in the Sirente-Velino Regional Park: Circuitscape and Omniscape implementations.

Here, I assess the overall landscape connectivity for the Mariscan brown bear (Ursus arctos marsicanus) in the Sirente-Velino Regional Park (SVRP) in central Italy http://www.parcosirentevelino.it/. In particular, I derive the energy landscape for the area and use it as a “resistance matrix” for animal movement, i.e. assuming that the animal would move across low energy landscape trajectories in order to minimize energy costs.

The required input for enerscape using the ARC model are:

  • a digital elevation model (DEM) of the area of interest (m)
  • body mass of the animal (kg)

I load the DEM and the park polygon data:

dem <- raster("../data/dem-vignette.tif")
park <- st_read("../data/SVRP-polygon.shp")
## Reading layer `SVRP-polygon' from data source `/home/eb97ziwi/Proj/enerscape-paper/data/SVRP-polygon.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 30 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 857140.5 ymin: 4665004 xmax: 897975.7 ymax: 4691681
## CRS:            unknown
m <- mask(dem, park, inverse = TRUE)
plot_area(dem, park, m, void = TRUE)
title("Elevation (m)")
scalebar(10000, 
         type = "bar", 
         divs = 2, 
         label = as.character(c(0, 5, 10)),
         below = "kilometers")

I calculate the energy landscape for a bear of 140 kg:

en <- enerscape(dem, 140, "kcal")
##  - Raster cells are assumed to have same horizontal and vertical resolution and with planar coordinate reference system (e.g. UTM)
##   | Calculating slope
##   | Calculating work
##   | Calculating conductance (1 / work)
##  - Do not use slope with negative values for distance calculations
plot(en$rasters, axes = FALSE, box = FALSE)

I choose two points for which to compute the least-cost path and the connectivity using https://docs.circuitscape.org/Circuitscape.jl/latest/usage/:

p <- data.frame(x = c(877367, 882653),
                y = c(4674192, 4677413))
plot_area(dem, park, m, void = TRUE)
points(p, pch = 20)

lcp <- en_lcp(en, p[1, ], p[2, ])

circuitscape_skeleton(en,
                      file.path("/home",
                                "eb97ziwi",
                                "Proj",
                                "enerscape-paper",
                                "data"),
                      p)

And I run Circuitscape in Julia.

julia> using Omniscape
julia> run_omniscape("/home/eb97ziwi/Proj/enerscape-paper/data/circuitscape.ini")
## [ Info: 2020-12-19 08:31:24 : Precision used: Double
## [ Info: 2020-12-19 08:31:27 : Reading maps
## [ Info: 2020-12-19 08:31:28 : Resistance/Conductance map has 186803 nodes
## [ Info: 2020-12-19 08:31:36 : Solver used: AMG accelerated by CG
## [ Info: 2020-12-19 08:31:36 : Graph has 186803 nodes, 2 focal points and 1 connected components
## [ Info: 2020-12-19 08:31:36 : Total number of pair solves = 1
## [ Info: 2020-12-19 08:31:37 : Time taken to construct preconditioner = 1.683233535 seconds
## [ Info: 2020-12-19 08:31:37 : Time taken to construct local nodemap = 0.008903604 seconds
## [ Info: 2020-12-19 08:31:37 : Solving pair 1 of 1
## [ Info: 2020-12-19 08:31:38 : Time taken to solve linear system = 0.526888286 seconds
## [ Info: 2020-12-19 08:31:39 : Time taken to write current maps = 0.559248897 seconds
## [ Info: 2020-12-19 08:31:39 : Time taken to complete job = 15.111858446
## 3×3 Array{Float64,2}:
## 0.0   1.0      2.0
## 1.0   0.0     31.8625
## 2.0  31.8625   0.0
cs <- raster("../data/circuitscape_cum_curmap.asc")
cs <- log10(cs)
cs <- cs - min(values(cs), na.rm = TRUE)
cs[cs >= quantile(cs, 0.999)] <- quantile(cs, 0.999) #remove far outliers
cs <- cs / max(values(cs), na.rm = TRUE)
plot_area(cs, 
          park, 
          m, 
          col = viridis::inferno(100), 
          void = TRUE)
points(p, pch = 3)
lines(lcp$Paths)
title("Connectivity between locations")

To compute the overall landscape connectivity, I use the https://docs.circuitscape.org/Omniscape.jl/stable/ algorithm, which runs Circuitscape iteratively on a moving window. I initialize the omniscape.ini file to be passed to Julia:

omniscape_skeleton(en,
                   file.path("/home",
                             "eb97ziwi",
                             "Proj",
                             "enerscape-paper",
                             "data"),
                   radius = 10)

And I run the Omniscape algorithm in a Julia terminal. Don’t forget to specify the number of parallel threads: https://docs.circuitscape.org/Omniscape.jl/stable/usage/#Parallel-Processing.

julia> using Omniscape
julia> run_omniscape("/home/eb97ziwi/Proj/enerscape-paper/data/omniscape.ini")
## Starting up Omniscape. Using 7 workers in parallel. Using double precision...
## Solving moving window targets...
## Progress: 100% | Time: 0:03:14
## Done!
## Time taken to complete job: 196.3408 seconds
## Outputs written to
## /home/eb97ziwi//home/eb97ziwi/Proj/enerscape-paper/data/omniscape_1
## (Union{Missing, Float64}[0.5152914765395912 0.5884433266779532 … 0.719651357069539 0.8600629588479131;
## 0.4834003129933446 0.6594116484239914 … 0.635221929581483 0.7110678186022065;
## … ;
## 0.2127925084620419 0.2067469779201847 … 0.39577453680538544 0.46681040575465904;
## 0.17382917047741175 0.2093726458198999 … 0.31763823508216504 0.3738227777731409]
## Union{Missing, Float64}[1.0885252650071 1.136985249629082 … 1.169180462521136 1.1104937306261533;
## 1.053217843455199 1.0090282247541849 … 0.9301240063409788 1.1859951520336032;
## … ;
## 1.1798282079383584 0.9871029439774996 … 0.9873244151514602 1.2432973368603717;
## 1.1120398150491893 1.168088759875787 … 1.0583987239714665 1.0877633126993709])

Sometimes, few cells with extremely high current/connectivity can make the output map hard to read, as most of the connectivity values are squeezed into a narrow range. To avoid this, the value of extremely high outliers can be forced to an upper limit. For example, I forced outliers to the 0.999 quantile value of the connectivity and standardized the omniscape output between 0 and 1.

os <- raster("../data/omniscape_1/cum_currmap.tif")
os <- os - min(values(os), na.rm = TRUE)
os[os >= quantile(os, 0.999)] <- quantile(os, 0.999) #remove far outliers
os <- os / max(values(os), na.rm = TRUE)
plot_area(os, 
          park, 
          m, 
          col = viridis::inferno(100), 
          void = TRUE)
title("Landscape connectivity")

References

Etten, Jacob van. 2017. “R Package Gdistance: Distances and Routes on Geographical Grids.” Journal of Statistical Software 76 (1): 1–21. https://doi.org/10.18637/jss.v076.i13.

Pontzer, Herman. 2016. “A Unified Theory for the Energy Cost of Legged Locomotion.” Biology Letters 12 (2): 20150935. https://doi.org/10.1098/rsbl.2015.0935.

Prampero, P. E. di, G. Cortili, P. Mognoni, and F. Saibene. 1979. “Equation of Motion of a Cyclist.” Journal of Applied Physiology 47 (1): 201–6. https://doi.org/10.1152/jappl.1979.47.1.201.